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Abstract: A reduced-order model algorithm, based on approximations of Lax pairs, is proposed 
to solve nonlinear evolution partial differential equations. Contrary to other reduced-order meth- 
ods, like Proper Orthogonal Decomposition, the space where the solution is searched for evolves 
according to a dynamics specific to the problem. It is therefore well-suited to solving problems with 
progressive waves or front propagation. Numerical examples are shown for the KdV and FKPP 
(nonlinear reaction diffusion) equations, in one and two dimensions. 
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Reduction de modeles basee sur des paires de Lax 

approchees 

Resume : Nous proposons une methode de reduction de modeles, basee sur des paires de 
Lax approchees, pour resoudre des equations devolution non lineaires. Contrairement a d'autres 
methodes de reduction de modeles, comme la POD, l'espace dans lequel la solution est cherchee 
evolue selon une dynamique reliee au probleme. La methode est par consequent bien adaptee a 
des problemes comportant des ondes progressives ou des propagations de fronts. Nous montrons 
des exemples numeriques pour les equations de KdV et FKPP en dimension un et deux. 

Mots-cles : Reduction de modele, paires de Lax, transformee de scattering, KdV, FKPP 
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1 Introduction 

This work is devoted to a method for solving nonlinear Partial Differential Equations (PDEs) by 
a reduced-order model (ROM) based on the concept of Lax pairs. 

In his seminal work [10 , Lax proposed a formalism to integrate a class of nonlinear evolution 
PDEs. He introduced a pair of linear operators C[u) and A4(u), where u denotes the solution of 
the PDE. The eigenfunctions of C(u), are propagated by a linear PDE involving Ai{u). For some 
PDEs. the eigenvalues of C{u) have the remarkable property to be constant in time, as soon as 
C(u) and Ai(u) satisfy the Lax equation (see ^ below). This formalism, which is closely related 
to the inverse scattering method, can be applied to a wide range of problems arising in many 
fields of physics, such as Korteweg-de Vries (KdV), Camassa-Holm, Sine-Gordon or nonlinear 
Schrodinger equations. 

The main idea of the present paper is to use the eigenfunctions of C{u) as a basis to ap- 
proximate the solutions of the PDE. Contrary to standard ROM techniques, like the Proper 
Orthogonal Decomposition (POD, see e.g. [TS]) or the Reduced Basis Method (RBM, see e.g. 
|12II13| ). the basis is therefore time dependent. This property makes the method especially well- 
suited to problems featuring propagation phenomena, which are known to be difficult to tackle 
with POD. 

To define a model that is genuinely "reduced", the number of eigenfunctions used to ap- 
proximate the solution has to be small. This seems to be the case for various problems of 
practical interest. This has been recently pointed out by Laleg, Crepeau and Sorine who used 
the eigenfunctions of a Schrodinger operator, playing the role of C{u), to provide a parsimonious 
representation of the arterial blood pressure [5JGEIE]- Their signal processing technique, called 
SCSA (for Semi-Classical Signal Analysis), has been the starting point of our work. 

In the huge literature devoted to Lax pairs (see for example |3j 0] and the reference therein) , 
operator M.(u) is generally used, or searched for, in closed- form. One of the contributions of 
the present work is to propose an approximation of M. (u) that can be used even when the Lax 
pair is not explicitely known. In addition, whereas most of the studies consider one dimensional 
domains and functions rapidly decreasing at infinity or periodic boundary conditions (see |2] for 
a theory on the finite line), our approximation method can easily be used on multi-dimensional 
bounded domains, with standard boundary conditions. 

The structure of the work is as follows: in Section [2] the Lax pair is introduced for the KdV 
equation and the link with the inverse Scattering Transform method is recalled; in Section [3] 
numerical approximations of the Lax operators are proposed; Section|4]is devoted to our reduced- 
order model algorithm, that will be called ALP, for Approximated Lax Pair. Some numerical tests 
are presented in Section [5] for the Korteweg-de Vries and Fisher- Kolmogorov-Petrovski-Piskunov 
equations in one and two dimensions. 



2 The inverse scattering method for nonlinear PDEs 

Lax method has strong connections with the inverse scattering method, which can be roughly 
viewed as the counterpart of the Fourier transform for nonlinear problems (see e.g. [TJ[3]). In 
this section, we briefly recall how the scattering transform is classically used to solve nonlinear 
PDEs. This reminder is not strictly necessary to introduce our method, but the reader might find 
it useful to understand the underlying ideas. We consider the standard example of Korteweg-de 
Vries equation (KdV): 

d t u + 6ud x u + d xxx u = 0, (1) 
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for (x,t) G IR x IR_|_, with the initial condition u(x,0) = uq(x),x G R. The Lax pair associated 
with this equation is defined by 



C(u)v = -d 2 xx v - uv, 
M(u)v — Ad xxx v + 6ud x v + 3vd x ii 



, .. , ( 2 ) 



These operators satisfies the Lax equation 

d t £(u)+£(u)M(u)-M(u)£(u) = 0, (3) 

if and only if u(x, t) is solution to the KdV equation. Consider the spectral problems, parametrized 
by t, 

C(u(x,t))4>(x,t) = X(t)4>(x,t). 
It can easily be deduced from the Lax equation that 

d t X = 0. (4) 

When this property is satisfied, the problem is said to be isospectral. The Lax equation describes 
how the operator £(u) evolves in time. Its eigenfunctions can be used to reconstruct the solution 
u at every time. 

The solution of the KdV equation by inverse scattering can be decomposed in three steps. 
The first step is to solve the spectral problem associated with the initial data: 

C{uq)<P = \<j) 

which is the linear Schrodinger equation with potential —uq. In general, this problem, known 
as the Direct Scattering Transform, has a continuous spectrum of positive values and a discrete 
spectrum of negative eigenvalues. For the sake of simplicity, we assume that uq is such that it can 
be reconstructed with the sole discrete part of the spectrum (such a uq is known as a reflectionless 
potential). We also assume that uq is in the Schwartz space <S(R), i.e. decays sufficiently rapidly 
at infinity, which implies that the number of negative eigenvalues is finite. For m = 1, • • • , iV_, 
we denote by <fi m the eigenfunctions normalized in L 2 (R), by A m the eigenvalues, we define 
Km = \/— A m and we rank the eigenvalues such that k\ > K2 > ■ ■ ■ > kjv_ ■ The Deift-Trubowitz 
formula then gives (see [I]): 

u o(x) = ^2 K ™4> 2 m {x,Q)- (5) 

m— 1 

When x goes to infinity, the eigenfunctions 4> m (x, 0) can be proved to be equivalent to c m (0) exp(— n m x). 
The quantities (k to (0), c m (0)) are known as the scattering data. Note that in general, the scat- 
tering data also contain the so-called reflexion coefficient, which is not present here because of 
the assumption made on uq. This missing part is important for scattering problems, but not for 
the present work. 

The second step consists in propagating the scattering data in time. This is trivial for the 
K m , since they are already known to be independent of t (see Q). For c m , it can be shown, 
using the fact that the solution is in the Schwartz space, that 

dr 

-^(t)=4K 3 mCm (t), (6) 

which readily gives c m (t). 



Inria 



Reduced-Order Modeling based on Approximated Lax Pairs 



5 



The third step is the Inverse Scattering Transform, i.e. compute u(x,t) from the scattering 
data (/t TO , c m (t)). In general, this can be done by solving an integral equation, known as the 
Gelfand-Levitan-Marchenko equation. When uq is reflectionless, i.e. given by §5§, a closed-form 
expression of u{x, t) using the scattering data is known. More details can be found e.g. in [TJ [3] . 

The method that has just been described heavily relies on the specificity of the KdV equation 
set in the whole space R, with an initial condition rapidly decreasing at infinity. In the remainder 
of this article, we propose a method directly based on an approximation of the Lax pair. It 
formally follows the same steps as the inverse scattering method, but can be used in a more 
general setting. 

3 Preliminary results 

We consider an evolution PDE in a bounded subset Q, of R d , d > 1, of the form: 

d t u = F(u), (7) 

where F(u) is an expression involving u and its derivatives with respect to x\,...,x&. The 
problem is completed with an initial condition 

u(x, 0) = Uq(x), for x = (xi, . . . , Xd) € fi- (8) 

For the sake of simplicity in this Section, u is assumed to vanish on the boundary dfl. Other 
boundary conditions will be considered in the numerical tests. 

The solution to this problem is searched for in an Hilbert space V (of functions vanishing on 
dfl) and approximated in Vh, a finite dimensional subspace of V, for example obtained by the 
finite element method. Let (vj h)j=i..N h denotes a basis of Vh and (•, •) the L 2 (f2) scalar product. 

We follow the three steps presented in the previous section to solve the KdV equation by 
inverse scattering. For each step, we propose a numerical approximation that will be used for 
the reduced order model integration. 



3.1 Semi Classical Signal Analysis 

Consider a nonnegative signaj^i(a;), for x — (x%, . . . ,Xd) € 51, a real number x > 0, and the 
linear Schrodinger operator 

C x {u)<f> = -A<t>-xu<f>, (9) 

where A denotes the Laplacian in d dimensions. The approximated scattering transform we 
consider is called SCSA for Semi Classical Signal Analysis. It has been proposed in [B], analyzed 
in [5] , and successfully used for different applications to signal analysis in hemodynamics [S] • 
It partially relies on the results proved in IH] and it consists of only keeping the modes 
corresponding to the negative eigenvalues (A„)„ = i.. jv_ to approximate u by the Deift-Trubowitz 
formula: 

JV_ 

U(x) = X^ 1 ^2 Km< ^m> ( 10 ) 

with K m = y'-Am. The parameter x > is chosen in order to reach the desired accuracy. For 
large values of x > 0, the representation is more accurate, but also more expensive since the 
number of negative eigenvalues is larger. 

If u(x) is not nonnegative, it is replaced by u(x) — min^gn uq(x) 
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Numerically, we search for <j) m ,h € Vh and X mt h such that 

<V<?!» m ,fc, Vv it h) - x(u4>m,h,Vi t h) = Ki,h{(j> m ,h,v l . h ) 1 for i= 1, .. . , iV/j, (11) 
and it is approximated by 

N- 

where \ IS chosen such that 

\\u - %|U=(n) < e , 

where e Q is a prescribed tolerance. 



Remark. Note that problem (111 implies: 



(u4> m .h,4>p,h) - -{^4>m,h^4>p.h) + — (^/i^p.h) = 0) for ^ = 1, ... , JV_ , 

X ' X 

which is the Euler-Lagrange equation of the maximization of the L 2 projection of u onto the 
space spanned by the § 2 m ft , up to a regularization term, and under the constraint of a unitary 
L 2 norm. 

3.2 Modal approximation of the Lax propagation operator 

For an arbitrary PDE Q, the construction of C(u) and M.(u) satisfying 



C{u)4> = A0, 
d t 4> = M(u)q 



(12) 



and the Lax equation |3| is not obvious. In addition, the isospectral property Q, which was 
instrumental in solving the KdV equation, is in general not satisfied. 

In this work, we choose C{u) as the linear Schrodinger operator associated to the potential 
— \u, as in ([9]). The following proposition shows that it is possible to compute an approximation 
of A4(u) in the space defined by the eigenfunctions of C(u) and to derive an evolution equation 
satisfied by the eigenvalues of C(u), even when the operator M(u) is not known in a closed-form. 

From now on, the eigenfunctions of C(u(t)) will be denoted by (</> m (i))m=i---A/_ when they are 
used to approximate the solution u(t), and by (ip m (t))m=i---N M when they are used to approx- 
imate the operator Ai(u(t)) or the evolution equation satisfied by the eigenvalues. Functions 
<j) m and i\) m are therefore the same, but the numbers iV_ and Nm will be different in general: 
while the solution will be approximated only on modes corresponding to negative eigenvalues, 
the approximation of the operators will also necessitate some modes associated with positive 
eigenvalues. Thus Nm will be larger than iV_ in practice. 

Proposition 1. Let u be a solution of equation Let C x (u) be defined by 

C x (u)ifj = — Atp — xuip (13) 

where x is a given positive real number. 

Let Nm G IN*. For m e {1, . . . , Nm}, let X m (t) be an eigenvalue of C x (u(x,t)), and ipm(x,t) 
an associated eigenfunction, normalized in L 2 (fl). Assume there exists an operator A4(u) such 
that 

d t ip m = M(u)ip m . (14) 
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Then the evolution of X m is governed by 

d t X m = -X( F ( u )^m,1pm}, (15) 

and the evolution of ip m satifies, for p £ {1, . . . ,Nm}, 

(d t ip m ,ipp) = M mp (u), (16) 

with 

M mp (u) = —(F(u)tp m ,ip P ), ifp^mandX p ^X m , 

Ap — A m (17) 

M mp (u) = 0, if p = m or X p = X m . 

We will denote by M(u) £ ^ n m*n m ^ g s k ew _ S y mme t r i c matrix whose entries are defined by 
M mp (u). 

Proof. Differentiating with respect to t the equation satisfied by the m-th mode 

C x (u(x,t))ip m (x,t) = \ m (t)i/) m (x,t), 

we get 

(C x (u) - X m T) dtipm = d t X m i>m + xF(u)ip m . (18) 
The scalar product is taken with a generic ip p , leading to: 

((C x (u) - X m X) d t ip m ,4> P } = d t Xm(^miipp) + x(F(u)ip m ,4> p ), (19) 

Using the self-adjointness of the operator and the orthonormality of the eigenfunctions, the 
following problem is obtained: 

(X p - X m ) (d t if> m , 4>p) = d t X m S mp + x(F(u)ip m ,ip p ). (20) 



Taking p — m, this proves (15 1. In addition, the L? norm of ip m being 1, (dtip m ,ipp) = 0, i.e. 

©2- 

If p =/= m, but X p — X m (multiple eigenvalues), we arbitrarilly set M mp (u) = 0. 
For X p ^ X m , we deduce from p0|: 



(dt^m^ P ) = X , (F(u)^ m ,i> p ), (21) 

Ap 



Combining (14 1 with (211, equation (17 1 is obtained, which completes the proof. <) 

Equation (17) gives an approximation of the operator M(u) on the basis defined by the 
modes at time t. This representation is convenient from a computational standpoint since it 
can easily be obtained from the expression F(u) defining the PDE Q, without any a priori 
knowledge of A4(u). With this approximation of Ai(u), the evolution of the modes can be 
computed according to the nonlinear dynamics of the system. This is an important difference 
with standard reduced-order methods, like POD, where the modes are fixed once for all. 

To set up a reduced order integration method, only a small number Nm of modes will be 
retained. This number has to be chosen in order to represent the dynamics in a satisfactory way. 
A possible indicator of the quality of the approximation is given by the following quantity 

N M 

e(ip m (t),N M ) = ^(Af m „(u(t))) 2 , (22) 

n=l 
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which is an approximation of the L 2 norm of the time derivative of ip m : 

„ N M N m 

/ (<9 f -0m) 2 dQ sa ^2 M m i(u)M mn (u)(ip n ,ipi) = ^ M mn {uf = e(ip m ,N M )- 

•* n n.l=l n=l 

By summing up over the modes, the Frobenius norm of the representation of the evolution 
operator is recovered: 

N M N u 

e(<M*), N M ) = (M mn (u(t))) 2 = F 2 Nm . (23) 

m m.n—1 

This norm may be used as an error indicator for the dynamics recovery and it was investigated 
by the numerical experiments presented in Section [5] 

3.3 Approximated Inverse Scattering Transform 

Proposition [T] gives an approximated way to propagate the eigenmodes and the eigenvalues 
associated to a Lax pair. From a given set of eigenmodes and eigenvalues, we will have to 
reconstruct the solution. We therefore need an approximated counterpart of the inverse scattering 
step presented in Section [2] 

Assume that a family (<j) m )m=i..N- of eigenfunctions of C x (u), and eigenvalues (A m ) m= i..jv_ 
are known. We propose to approximate u as 

N— 

u = 

k=l 

Inserting this expression in 

(c x (u)4> m , 4> p ) = \ m {<t> m , <t>p) , 
and using that (<f> m , 4> P ) — 5 mpi we obtain: 

N - ! 

= (\ m + {A<!>m,<t>m})- (24) 

fe=l A 

This is a linear system for ctk, whose resolution is costless when a small number of modes is 
considered. 



4 Reduced- Order Modeling based on Approximated Lax 
Pairs (ALP) 

The eigenvalues and the eigenfunctions of C x (u(-,t)) evolve as t changes. We propose in this 
section numerical schemes to approximate these evolutions. 

4.1 Approximated propagation of the eigenvalues 

According to Proposition [l] the evolution of the eigenvalues is governed by (15 1. Numerically, 
we consider a simple explicit Euler scheme: 

Assume that u n , A™ and (ipp)p=i,,,.,N M , are known at step n, compute A™ +1 as: 

x; +1 = x;~5t(F{u n )r p ,r P ), (25) 

where St denotes the time step. 
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4.2 Approximated propagation of the eigenfunctions 



To propagate the modes ip p , it is important to devise a procedure that preserves their orthonor- 
mality. We first note that equation (16 1 yields: 



dt 



N M 

= ^M mp (u)il) p 

p=l 



(26) 



where r m (t) e [span (V>i (t), . . .,ip NM (t))} . 

Assume that u n and (ipp) p =i t ...,N M are known and that ("0m> V'p) = ^mp for m,p — 1, . . . , Nm . 
Consider the following problem, obtained from (26 1 by neglecting the residual r m and linearizing 
the first term in the right-hand side: 



dlpr, 



P =i 



whose solution is 



Note that 



dt 

N M 

^ m (t) = ^[exp(M(u")t)] mp V p l . 
P =i 

$ m (t),Mt)) = [exp(M(u n )t)exp(M(u n ) T t)] 



m p 



(27) 



(28) 



(29) 



since M{u n ) is skew-symmetric. Thus, by construction, this procedure preserves the orthonor- 
mality. Ideally, V>m +1 could be defined as ip(6t). But to avoid the expensive computation of the 
exponential of a matrix, we suggest to approximate exp(M(u n )5t) by T(u n ) = I + M(u n )5t + 
\M(u n ) 2 8t 2 . Thus, is eventually defined by: 



N M 



St 2 

I + M (u n )5t + —M(u n ) 2 



With this approximation, the orthonormality of ('0m +1 )m=i....,Ar J 
error of the order of St 4 . 



(30) 

is only perturbed with an 



4.3 The ALP algorithm 

We now have all the tools necessary to the definition of a reduced order model based on the 
approximated Lax pairs. 

Initialization: Let uq be the initial condition and let eo > be a prescribed tolerance. 
Compute a set of modes (ij;^ h ) m =i...N M and eigenvalues (A^ h ) m =i...N M by solving 

(W?™^, Vv i>h ) - X{uipm,h: v i,h) = >>m,h('>Pm,h> v i,h), for « = 1, . . .,N h , 

The eigenvalues A^ h are ranked in the ascending order. Let < Nm be the number of negative 
eigenvalues. Throughout the algorithm, we define (j) 1 ^ h = ip" h , for m G {1, . . . , iV_}. The initial 



condition is approximated by the SCSA method (see Section 3.1) 
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where x is chosen such that ||uq — u^JIl 2 ^) < eo. 

Time evolution: For n > 0, assume that u 1 ^, (ip^ h ) 



-i-..jsr M and Aj! are known at time t n . 



1. Compute the matrix M{u™ L ) approximating the operator Ai(u(t n )) with (17) 
M mp «) = ^(^(OC^^J- for m,p = 1,.. .,N M . 



2. Compute the eigenvalues with (|25|) 



3. Compute the eigenfunctions tp^ 1 with (30 1: 



N M 

<l = E 

P =i 



St 2 

I + 6tM(ut) + —M(u n h ) 2 



fl>. 



4. Solve system ([24} for a™ +1 



JV_ 



E< +1 (^m>™,^) = -\ ( a p +1 + 



P =i 



5. Compute u? +1 with 



p=l 



r 1 (^i 1 )' 



Remark. For the sake of simplicity, N- and Nm are fixed in this algorithm. They of course 
could be adapted along the resolution acording to various criteria. We will not investigate this 
possibility in the present study, except in the following case: if one eigenvalue which was initially 
positive becomes negative, then the corresponding eigenmode can be simply added to the set 
used to approximate the solution, and the value of 7V_ incremented accordingly. This will be 



used in the test case presented in Section 5.2.1 



5 Numerical Experiments 

In this section, some numerical experiments are presented. The aim is to assess the proposed 
technique and to highlight some differences with the POD. 

Two test cases deal with the Korteweg-de Vries (KdV) equation. As recalled in Section [2| 
this equation is a remarkable example of an integrable system, which can be exactly solved by 
inverse scattering (see OH]). This problem is therefore a good candidate to assess our numerical 
approach in a case when an analytical solution is possible. 

Then, the Fisher- Kolmogorov-Petrovski-Piskunov (FKPP) equation is considered in one and 
two dimensions. This equation, which arises in many applications, features front propagation. 
Contrary to KdV, it is not isospectral, and no closed-form expression of the Lax pair is known. 
So, it is an interesting problem to test our method when no analytical approach is available. 
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(a) 



(b) 



Figure 1: a) Exact one-soliton solution, see Eq.(31 ), at initial and final time (T=5.0) b) Relative 
error for the Frobenius norm of the operator M and for the peak position, at t = T, in semi- 
logarithmic plot. 



5.1 Korteweg-de Vries equation 

We consider the one-dimensional KdV equation , with an initial condition uq that will be either 
a one-soliton or a three-soliton. The computational domain is bounded for practical reason, but 
the time range of the computations and the length of the domain are chosen such that the 
boundaries have no effect on the solution. 



5.1.1 One-soliton propagation 

The reference solution shown in Fig(l]a) at initial and final time (T = 5.0) is a one-soliton 
propagation, namely: 

u(x,t) = ^sech 2 (^(x-0t)j, (31) 

with p = 4. 

The modes are extracted by using the initial condition only uq = u(x, 0), setting x = 1. The 
Schrodinger spectral problem was discretized in a space of Nh — 500 piecewise linear functions. 

As uq is the initial datum of the one-soliton propagation and \ = 1 provides the analytical 
expression for C(u) in the case of the KdV equation (see Eq.Q), only one eigenvalue belongs 
to the discrete spectrum and the corresponding mode squared is exactly uq, up to discretization 
errors (10 -4 in L? norm for the present case). Hence, only one mode is retained in order to 
represent the solution (iV_ = 1). 

The influence of the number Nm of modes on the dynamics approximation is investigated by 
computing the error in L? norm of it as a function of time. The results are shown in Fig(3]a), 
in a semi-logarithmic plot. As expected, the error globally decreases when the number of modes 
is increased and it tends to increase in time, due to the accumulation of the errors during the 
integration. 

A meaningful and synthetic criterion of the dynamics recovery quality is the relative error 
in the peak position at final time, that is, the difference, in modulus, between the position of 
the peak of the exact and the reconstructed solution, normalized by the distance traveled by the 
wave, which is 20 in this case, see Fig. Ilia). 
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(a) 



(b) 



Figure 2: Comparison between the exact solution and the solution obtained with ALP when a) 
10 modes and b) 25 modes are used in order to represent the evolution operator A4(u). Dashed 
line represents the initial condition. 



It is enlightening to compare the error in the peak position and an indicator of the approx- 
imation error of operator AA(u). To estimate this indicator, we compute the relative differ- 
ence between a reference value and the Frobenius norms (see Section 3.2) of matrices M(u) for 
N M = 10, 15, 20, 25, 50, computed at final time t = T: 



£ F (N M ,N a 



F, 



F 



(32) 



where Fn — (J2 m p=i ^mp) 1 ^ 2 an d where the reference value Nqq = 100 (adding more modes 
changes its value less than 0.5%). In Fig. [l]b), the error in the peak position and £f(Nm, ^oo) 
are shown as a function of Nm, in & semi-logarithmic plot. The trend of the curves is practically 
the same, suggesting that £f(Nm, -Woo) might be used as an indicator to assess the quality the 
approximation. 

As said above, the KdV equation defines an isospectral flow, that is, the eigenvalues of the 
Schrodinger operator associated with the solution of the problem do not vary in time. An 
interesting validation test is to see how far the reduced order approach satisfies the isospectrality 



property. The computation of the right-hand side of Eq.( 15 I provides the right result, for all the 
number of modes used. When St — 2.5 10~ 3 is used for the time discretization, the difference 
between the right value of the first (negative) eigenvalue (namely A = — 1) and the computed 
one is 4.5 10" 5 at final time t — T. 

Let us recall that only one mode is used in order to represent the quantity u, while a varying 
number of modes is adopted to study the dynamics. In Fig|2]a comparison between the exact 
(see Eq.(31l) and the reduced-order integration method solutions is shown. In Fig{2ja) the final 
time solution is compared to the one obtained by representing the evolution operator M (u) with 
10 modes: the peak position, the amplitude and the wave shape are not correctly represented. 
Instead, when 25 modes are used, all the feature of the solution are well rendered, as it may be 
seen in Fig|2]b). 



In Fig|3lb) the portrait of the operator Ai(u) at final time is shown. According to (16), the 
entry M mp {u) represents the projection of the time derivative of the to— th mode onto the p— th 
one, at the current time. In this case, the highest values are all concentrated on the projection of 
the time evolution of the first mode (the only one used to represent u) on the others. Furthermore, 
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(a) (b) 

Figure 3: a) Errors in L 2 norm while varying the modes number in the approximation of the 
operator A4(u), b) representation of A4(u), in absolute value. 



the entries become smaller and smaller as m,p increase, which is desirable when setting up a 
reduced order integration approach. 

Remark. In this presentation, we restricted ourselves to PI finite elements for the sake of con- 
ciseness. During our investigations, we also used P2 finite elements and Hermite polynomials, 
and we noticed that much less modes were necessary to achieve a good approximation in those 
cases. The link between the quality of representation by the modes and the quality of approxima- 
tion of the underlying discretization method is an interesting question that might be addressed 
in future works. 



5.1.2 Three-soliton propagation 

In this section, we consider the propagation of a three-soliton by the KdV equation. The reference 
solution, shown in Figj4]a-b), has been generated by considering the Gelfand-Levitan-Marchenko 
equation, that, when solved for the KdV equation (see [T]) provides: 

u(x, t) = -2d 2 x log(det(7 + A(x, t))), (33) 
where A 6 R nx ™ [ s the interacting matrix, written in terms of the scattering data. In particular: 
A mn (x,t) = ° m 2 n exp{(fe m + k n )x - (k 3 m + kl)t) , (34) 

where fc m , c„ are 2n scalar parameters that may be linked to position and speed of solitons (see 
HHS]). For the present case: c= [5.0 1CT 2 , 1.5 10-\l.O 10 1 ], k = [1.0,1.5,1.75], i€ (-15,15) 
and t £ (0,0.5). An interesting feature of this test case is the complex interaction between the 
waves: at the end of the simulation, two of them are fused together (see Figjljb)). 

The scattering transform is applied at initial time and, by setting \ = 1, three distinct 
negative eigenvalues are found. This is in agreement with the analytical results and highlight 
the ability to decompose a traveling (nonlinearly interacting) waves system in modes which 
propagates separately. Thus, three modes are retained to represent the solution. 

To represent the evolution operator A4(u), we take Nm = 20 modes. With this value of Nm, 
we checked at every time step that the relative difference in Frobenius norm £f{Nm — 1,Nm) 



(see (32 1) was less than 1.0%. For the present case, the entries in the first three rows of its 
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(a) 



(b) 



(c) 



Figure 4: Plot of the exact three-soliton solution at a) initial and b) final times. In c), comparison 
between exact solution and reduced order model, at final time, when 20 modes are used to 



represent the operator. Problem setting defined by Eq.(33l 





(a) 



(b) 



Figure 5: Comparison between POD and the ALP a) initial time, static reconstruction with 3 
modes for both the techniques b) final time, integrated result with ALP, best possible recon- 
struction using 20 POD modes. 



representation, corresponding to the projection of the time variation of the three waves used to 
recover u on all the others, are higher in absolute value. It has been observed in the numerical 
experiments that, when k modes are necessary and sufficient to have a good representation of 
u(x, t), the matrix M has the first k rows (and columns, M being skew-symmetric) which contains 
the entries that contribute, maximally, to the Frobenius norm of the operator. 

In Figgb) a comparison at final time between the exact (see Eq.(|33[)) and the ROM solutions 
is shown. The dynamics, in particular the coalescence phenomenon, is correctly represented. 

Compared to the POD method (Proper Orthogonal Decomposition, see e.g. [E]), an inter- 
esting feature of the ALP is its ability to extrapolate traveling wave type solutions out of the 
database used to defined the approximation space. To illustrate this point, we compared the two 
methods on the three-soliton test case. A POD basis was built by using the Sirovich technique 
(see [TS] for details), by taking N s — 50 snapshots on half the time history, i.e. between t = 
andt = 0.25. InFigj5ja), we compare the ability of the two approaches to approximate the initial 
data with only 3 modes. In Fig|5]b) we compare the solution at final time t — 0.5. The POD 
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(a) (b) (c) 

Figure 6: Contour of the modes squared used to represent u for a three-soliton propagation (see 



Eq.(33l): a) first b) second c) third mode, 30 levels between maximum and minimum. 




(a) 



(b) 



(c) 



Figure 7: Contour of the modes for a three-soliton propagation (see Eq.(33)): a) fourth b) fith 
c) sixth mode. 30 levels between maximum and minimum. 



reduced order model was run with 20 modes, and the ALP with N— = 3 modes to represent the 
solution and Nm = 20 for A4(u). In this example, the accuracy of the ALP clearly outperforms 
the POD. Let us emphasize that POD was built by using half of the time history while the ALP 
simply propagates the initial datum, without any pre-computed database. 

In Fig(6]the contour of the modes squared used to represent the solution - i.e. the <j)f n used 
in (10 1 - is shown in t, x plane. There are three distinct waves, propagated separately. In Fig (7] 



the contour of the first three modes tp^ associated with positive eigenvalues is shown. These 
modes are not used to reconstruct the solution, but they contribute to the waves evolution, in 
the sense that the projection of the time derivative of the waves on this modes is significant. 

Remark. For the KdV equation, as well as for any equation admitting a Lax pair, the expression 
of the operator M. (u) is known analytically. This may be profitably used for the reduced order 
model, since in that case, there is no need to search for an approximated representation of the 
propagation operator. In other words, Steps 1 and 3 of the ALP algorithm can be replaced by 



the direct solution of equation (12)2- We checked for the KdV equation that the two approaches 
give similar results. 
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5.2 Fisher-Kolmogorov-Petrovski-Piskunov equation 

In this section the FKPP equation is considered as an example of non-isospectral flow equation, 
in a finite domain, with homogeneous Dirichlet or Neumann boundary conditions. This differs 
substantially from the setting usually adopted when the Lax pair is known in a closed-form. 

5.2.1 ID FKPP with homogeneous Dirichlet boundary conditions 

The equation reads: 

dtu = d^ x u + au(l — u), in ft = [0, 1] 
u = 0, on oil 

For the present case a = 1.0 • 10 3 , the space domain was [0,1] and Nf = 250 piecewise linear 
functions were considered. The time domain was [0, 7.5 • 10~ 3 ] and 100 integration points were 
taken. The reference solution was obtained by discretizing in space by means of piecewise linear 
functions and by using a mixed implicit-explicit scheme in time: the linear diffusion part of the 
equation was discretized by means of a Cranck-Nicolson scheme, the nonlinear term by an explicit 
second order Adams-Bashforth scheme, with St = 7.5 10~ 5 . In FigjSJa) the reference solution 
is plotted at different times: the initial solution is taken as: uo = exp (— 10 2 (x — 0.25) 2 ) + 
exp(-100(>-0.75) 2 ). 

The number of modes to be retained in order to represent u and A4(u) have to be chosen. The 
first one is determined by setting \, f° r which no theoretical value is available, contrary to what 
happened for the KdV equation. The scattering parameter \ — 5 ■ 10 2 is chosen as explained 
in the initialization stage of the ALP algorithm, with eo = 10~ 3 . This corresponds to 4 distinct 
modes <pk associated to negative eigenvalues. The number of ipk modes for the approximation 
of M(u) is Nm = 10, chosen with the same criterion used in the Section 5.1.2| (Frobenius norm 
modify by less that 1% when adding a new mode). 

In Figj8]b) a comparison between the reference initial solution and its reconstruction is shown. 

In Fig(9ja) a comparison between the reference solution at t = T and the reconstructed one 
is shown. The error is mainly due to the appearence of wavy structures. However, the dynamics, 
wich is non-trivial (coalescence of structures) is correctly recovered. 
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Figure 9: a) Comparison between the exact solution of Eq.(35l and the reconstructed one at 
final time; b) relative error in L 2 norm with respect to time 



In Fig|9jb) the relative error in L 2 norm is shown as a function of time. It is smaller than 10%, 
and in this case it is not monotonically increasing in time. This behavior is due to the eigenvalues 
dynamics. Contrary to what happened with the KdV equations, the eigenvalues are evolving in 
time and some positive eigenvalues may become negative. At the very beginning, A_ = 4 modes 
(corresponding to the negative eigenvalues) are sufficient to represent the solution. Although it 
is possible to keep these modes throughout the simulation, we noticed that the approximation 
was improved if, during the simulation, we added the modes corresponding to the new negative 

" 3 and 



4.3 



In Fig 9 b) the error peaks at times t = 1.5 • 10 



eigenvalues, as indicated in Remark 

t = 2.5 • 10~ 3 corresponds to the instant at which A5 and A6 respectively change their sign, and 
the corresponding <fi 5 and fa are used to approximate the solution. 



5.2.2 2D FKPP with homogeneous Neumann boundary conditions 

The 2D FKPP equation reads: 

d t u = Au + au(l — it), in f2, d n u — on dfl, 



(36) 



where a — 10 3 for the presented test case and homogeneous Neumann boundary conditions are 
imposed on dft, SIC K 2 . A reference solution for this problem was obtained by integrating the 
equation by means of a PI finite element method. The geometry of the domain and the mesh are 



shown in Fig 10a). The boundary is a set of line segments of unitary and integer multiple of unity 
length; the mesh counts 9490 triangles. In order to advance the initial solution in time a hybrid 
Crank-Nicolson, Adams-Bashforth 2 scheme was implemented, to discretize the diffusive and the 
reaction terms respectively. The time step was chosen as St = 2.5 10~ 4 . The initial condition 
a) ) was obtained by making a Gaussian u g = exp(— 50((a; — 2.5) 2 + (y — 0.5) 2 ) evolve 
75 time steps. The reference solution is featured by a nonlinear front propagation. 



(see Fig 
for n = 



11 



In particular, this geometry makes the evolution of this system particularly challenging to be 
recovered. First, the front propagates upward, then, once arrived at the T bifurcation it starts 



spreading almost spherically (see Fig 11 b)), and finally, due to Neumann boundary condition, 
when it reaches the upper wall, it splits and two fronts start propagating horizontally with ±x 



direction (see Fig 11 c) ). The overall evolution took n = 125 iterations. 

The initial datum was recovered by means of N- = 6 modes, obtained by solving the scatter- 
ing problem with x = 40. In Fig(l0]b) the contour of the initial relative error field is represented. 
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(b) 



Figure 10: a) Geometry of the domain Q and related mesh for Eq.(36l discretization b) Contour 
for the relative error field on the initial datum approximation, 25 isoline between maximum and 



minimum. 



(a) 




Figure 11: Contour levels of the reference solution of Eq.(|36j) on half the domain, 30 levels 
between maximum (red) and minimum (orange) at a) t=0, b) t=T/2, c=T. 
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(a) 




(c) 



Figure 12: Contour levels of the solution obtained by ALP, visualized on half the domain, 30 
levels between maximum (red) and minimum (orange) at a) t=0, b) t=T/2, c=T. 




Figure 13: Contour levels of the most significant mode (fith one) on half the domain, 30 isolincs 
between maximum and minimum at a) t=0, b) t=T/2, c=T. 



It is featured by oscillations, especially located near the lower bound of the domain, and it is of 
overall 10% in L 2 norm. It is due to the fact that the laplacian operator is not suitable for this 
kind of equation in finite domain, and the use of another compact operator should improve the 
quality of the datum reconstruction. The choice of the regularization operator is essentially a 
guess even in an analytical approach to the scattering transform (see [Tj for a detailed discussion 
on this topic). This point will be the object of further works. 

In Fig ]12| the contour levels of the solution obtained by propagating the modes are shown, 
at the same time instants of those ones shown for the reference solution, using the same scale. 
Despite the fact that the initial datum reconstruction is affected by some error, all the elements 
of the dynamics, that is the front speed and the global movement are recovered in a satisfactory 
way. Some improvement is needed on the accuracy in the solution representation (i.e. the front 
shape). 

This represents another case in wich POD would not be able to extrapolate out of database 
(by construction) and tends to perform poorly in terms of reduction. 

The modes were propagated by means of the proposed algorithm with a time step of 6t = 
1.0 10~ 3 . 



In Fig 13 a-c) the contour of the most significant mode is represented, at three different times, 
on half the domain (it is symmetric with respect to the y axis). The mode travels and then splits, 
its dynamics being similar to the one of the solution, featured by large displacements. 
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6 Conclusions and perspectives 

We have proposed a new reduced-order model technique, called ALP, consisting of three stages. 
First, a set of modes is constructed by the Semi Classical Signal Analysis method, i.e. as a set 
of orthonormal eigenfunctions of a linear Schrodinger operator associated to the initial condition 
(approximated scattering transform). Second, the modes are propagated by an approximation of 
a Lax operator, and the eigenvalues are updated. Third, the solution is reconstructed by solving 
a problem that plays a role similar to the inverse scattering transform. 

This approach allows us to set up a reduced order discretization that seems to be efficient 
for those problems involving progressive wave or front propagation. Contrary to other reduced- 
order methods, like POD, the solution can be extrapolated out of the database of pre-computed 
solutions. The method has been successfully tested on the KdV and FKPP equations, in ID and 
2D. 

The application of ALP to other problems is currently under investigation, in particular to a 
set of Euler equations modeling a network of arteries and to cardiac electrophysiology problems. 

Many questions would deserve further investigations: the number of modes used to approxi- 
mate the solution or the propagation operator could be adapted along the resolution, for example 



based on the indicator ( 32 ) ; other operators than the Laplacian might used for the scattering step; 



other schemes could be used to solve (27); positive eigenvalues might be used to approximate 



the solution, etc. This will be the subject of future works. 
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